
************************************
*Numeric Export disabled
local TEX ""
local TEXON 0 //Change "0" to "1" to enable exporting numeric results to .txt files
if !`TEXON' local TEX "*"
************************************
************************************
*Figure Export disabled
local FIG ""
local FIGON 0 //Change "0" to "1" to enable exporting figures (RELEVANT FOLDERS MUST FIRST BE CREATED)
if !`FIGON' local FIG "*"
************************************

*Necessary packages
*ssc install cem

************************************
*Set working directory
************************************

*Load Gilens' data, recode variables
do "Enns_SociologicalScience_Reproduction\LoadGilensDataRecodeConfirmAccuracy.do"


*******************************
*******************************
*Generate Figure after CEM
*Figure 3 (b and d)
*Figure 3 a and b reproduced from Figure 2 
*******************************
*******************************
****************
**Generate 90/10 variables to be used in cem program
****************

**90/10
quietly: gen dif_90_10=.
quietly: recode dif_90_10 .=1 if (pred90_sw - pred10_sw)>.1
quietly: recode dif_90_10 .=0 if (pred10_sw - pred90_sw)>.1

label var dif_90_10 ">10% dif b/t 90th and 10th income percentiles"

label define diflbl9010 0 "0 10th income percentile more supportive" 1 "1 90th income percentile more supportive"
label values dif_90_10 diflbl9010


************************************************
************************************************
*Plot distribution of preferences after CEM
*Figure 3 b (90th)
************************************************
graph set window fontface "Times New Roman"

*CEM
*set seed so figures are reproducable
set seed 32920671
quietly: drop if dif_90_10==.

******
*90th income percentiles
******
quietly: cem pred90_sw if policyadopted!=., tr(dif_90_10) k2k
*Change proportion variables to percent
gen pred90_sw_perc = 100*pred90_sw 
gen pred10_sw_perc = 100*pred10_sw 

****
*Generate predicted values after CEM
logit policyadopted pred90_sw_lor  if  abs(pred90_sw - pred10_sw)>.1 & cem_matched==1
mat define b = get(_b)
gen yhat90=exp(b[1,2]+(b[1,1]*pred90_sw_lor))/(1+exp(b[1,2]+(b[1,1]*pred90_sw_lor))) if e(sample)

logit policyadopted pred90_sw_lor  if (pred90_sw - pred10_sw)>.1 & cem_matched==1
mat define b = get(_b)
gen yhat90_greater=exp(b[1,2]+(b[1,1]*pred90_sw_lor))/(1+exp(b[1,2]+(b[1,1]*pred90_sw_lor))) if e(sample)
mat drop b

logit policyadopted pred90_sw_lor  if (pred10_sw - pred90_sw)>.1 & cem_matched==1
mat define b = get(_b)
gen yhat90_lesser=exp(b[1,2]+(b[1,1]*pred90_sw_lor))/(1+exp(b[1,2]+(b[1,1]*pred90_sw_lor))) if e(sample)
mat drop b

*Figure
local bins=40
twoway ///
(histogram pred90_sw_perc if (pred10_sw_perc - pred90_sw_perc)>10 & cem_matched==1, freq fcolor(navy%80) ///
lcolor(black) lwidth(vvthin) bin(`bins') yaxis(1))  ///
(histogram pred90_sw_perc if (pred90_sw_perc - pred10_sw_perc)>10 & cem_matched==1, ///
freq fcolor(orange%60) lcolor(black) lwidth(vvthin) bin(`bins') yaxis(1) ///
ytitle("Frequency", size(medlarge) axis(1)) ylabel(, nogrid notick axis(1))) ///
(scatter yhat90 pred90_sw_perc, yaxis(2) msize(medsmall) mcolor(black) mlcolor(none) ///
ytitle("Predicted Probability of Policy Change", ///
size(mlarge) axis(2)) ylabel(0(.2)1, nogrid notick axis(2))) ///
(scatter yhat90_greater pred90_sw_perc, yaxis(2) msize(medsmall) mcolor(orange%80) mlcolor(none)) ///
(scatter yhat90_lesser pred90_sw_perc,  yaxis(2) msize(medsmall) mcolor(navy%90) mlcolor(none)),  ///
xtitle("Percent Favoring Policy Change (90th Income Percentile)", size(medlarge)) ///
ytitle("Predicted Probability of Policy Change", size(large)) /// 
graphregion(fcolor(white) lcolor(white)) ///
plotregion(lcolor(black) lwidth(medium)) ///
xlabel(10[10]90, nogrid) xscale(lwidth(none) r(5 95)) ///
yscale(lwidth(none) r(0 1)) legend(off) ///
ysize(3) xsize(4)  

`FIG'graph export Figures/histogram90_cem.eps, replace
`FIG'graph export Figures/histogram90_cem.pdf, as(pdf) replace

drop cem_strata cem_weights cem_matched yhat90 yhat90_greater yhat90_lesser


***********************
*10th income percentiles
*Figure 3 c
***********************
quietly: cem pred10_sw if policyadopted!=., tr(dif_90_10) k2k

*Generate predicted values after CEM
logit policyadopted pred10_sw_lor  if  abs(pred90_sw - pred10_sw)>.1 & cem_matched==1
mat define b = get(_b)
gen yhat10=exp(b[1,2]+(b[1,1]*pred10_sw_lor))/(1+exp(b[1,2]+(b[1,1]*pred10_sw_lor))) if e(sample)

logit policyadopted pred10_sw_lor  if (pred90_sw - pred10_sw)>.1 & cem_matched==1
mat define b = get(_b)
gen yhat10_lesser=exp(b[1,2]+(b[1,1]*pred10_sw_lor))/(1+exp(b[1,2]+(b[1,1]*pred10_sw_lor))) if e(sample)
mat drop b

logit policyadopted pred10_sw_lor  if (pred10_sw - pred90_sw)>.1 & cem_matched==1
mat define b = get(_b)
gen yhat10_greater=exp(b[1,2]+(b[1,1]*pred10_sw_lor))/(1+exp(b[1,2]+(b[1,1]*pred10_sw_lor))) if e(sample)
mat drop b

*Figure
local bins=40
twoway ///
(histogram pred10_sw_perc if (pred10_sw_perc - pred90_sw_perc)>10 & cem_matched==1, freq fcolor(navy%80) ///
lcolor(black) lwidth(vvthin) bin(`bins') yaxis(1))  ///
(histogram pred10_sw_perc if (pred90_sw_perc - pred10_sw_perc)>10 & cem_matched==1, ///
freq fcolor(orange%60) lcolor(black) lwidth(vvthin) bin(`bins') yaxis(1) ///
ytitle("Frequency", size(medlarge) axis(1)) ylabel(, nogrid notick axis(1))) ///
(scatter yhat10 pred10_sw_perc, yaxis(2) msize(medsmall) mcolor(black) mlcolor(none) ///
ytitle("Predicted Probability of Policy Change", ///
size(mlarge) axis(2)) ylabel(0(.2)1, nogrid notick axis(2))) ///
(scatter yhat10_greater pred10_sw_perc, yaxis(2) msize(medsmall) mcolor(navy%90) mlcolor(none)) ///
(scatter yhat10_lesser pred10_sw_perc,  yaxis(2) msize(medsmall) mcolor(orange%80) mlcolor(none)),  ///
xtitle("Percent Favoring Policy Change (10th Income Percentile)", size(medlarge)) ///
graphregion(fcolor(white) lcolor(white)) ///
plotregion(lcolor(black) lwidth(medium)) ///
xlabel(10[10]90, nogrid) xscale(lwidth(none) r(5 95)) ///
yscale(lwidth(none) r(0 1)) legend(off) ///
ysize(3) xsize(4)  

drop cem_strata cem_weights cem_matched yhat10 yhat10_lesser yhat10_greater

`FIG'graph export Figures/histogram10_cem.eps, replace
`FIG'graph export Figures/histogram10_cem.pdf, as(pdf) replace

******************************************************************
*Generate Figures A-4: c (Using all data, not just data w/ preference gap)
*A-4 a from Figure 2
*A-4 b generated above
******************************************************************

************************************************
*10th (90 vs 10): Figure A-4 c
************************************************
*Need to re-generate predicted values using all the data (not just when preference gaps exist)
do "Enns_SociologicalScience_Reproduction\LoadGilensDataRecodeConfirmAccuracy.do"
*Change proportion variables to percent
gen pred90_sw_perc = 100*pred90_sw 
gen pred10_sw_perc = 100*pred10_sw 
logit policyadopted pred10_sw_lor  if  abs(pred90_sw - pred10_sw)>0
mat define b = get(_b)
gen yhat10=exp(b[1,2]+(b[1,1]*pred10_sw_lor))/(1+exp(b[1,2]+(b[1,1]*pred10_sw_lor))) if e(sample)
mat drop b

logit policyadopted pred10_sw_lor  if (pred90_sw - pred10_sw)>0
mat define b = get(_b)
gen yhat10_greater=exp(b[1,2]+(b[1,1]*pred10_sw_lor))/(1+exp(b[1,2]+(b[1,1]*pred10_sw_lor))) if e(sample)
mat drop b

logit policyadopted pred10_sw_lor  if (pred10_sw - pred90_sw)>0
mat define b = get(_b)
gen yhat10_lesser=exp(b[1,2]+(b[1,1]*pred10_sw_lor))/(1+exp(b[1,2]+(b[1,1]*pred10_sw_lor))) if e(sample)

local bins=40 
twoway ///
(histogram pred10_sw_perc if (pred10_sw_perc - pred90_sw_perc)>0, freq fcolor(navy%80) ///
lcolor(black) lwidth(vvthin) bin(`bins') yaxis(1)) /// 
(histogram pred10_sw_perc if (pred90_sw_perc - pred10_sw_perc)>0, /// 
freq fcolor(orange%60) lcolor(black) lwidth(vvthin) bin(`bins') yaxis(1) ///
ytitle("Frequency", size(medlarge) axis(1)) ylabel(, nogrid notick axis(1))) ///
(scatter yhat10 pred10_sw_perc, yaxis(2) msize(medsmall) mcolor(black) mlcolor(none) ///
ytitle("Predicted Probability of Policy Change", ///
size(mlarge) axis(2)) ylabel(0(.2)1, nogrid notick axis(2))) ///
(scatter yhat10_greater pred10_sw_perc, yaxis(2) msize(medsmall) mcolor(orange%80) mlcolor(none)) ///
(scatter yhat10_lesser pred10_sw_perc,  yaxis(2) msize(medsmall) mcolor(navy%90) mlcolor(none)),  ///
xtitle("Percent Favoring Policy Change (10th Income Percentile)", size(medlarge)) /// 
graphregion(fcolor(white) lcolor(white)) ///
plotregion(lcolor(black) lwidth(medium)) ///
xlabel(0[10]100, nogrid) xscale(lwidth(none) r(0 100)) ///
yscale(lwidth(none) r(0 1)) legend(off) ///
ysize(3) xsize(4) 

drop yhat10 yhat10_lesser yhat10_greater
`FIG'graph export Figures/histogram10_all.eps, replace
`FIG'graph export Figures/histogram10_all.pdf, as(pdf) replace


